Ejercicio propuesto por: Mariela Rajngewerc
La imagen original fue bajada de la página del earthexplorer. En esa página se pueden bajar imágenes con distinto nivel de pre-procesamiento. Para este ejercicio bajamos una imagen de nivel de procesamiento 2, esto quiere decir que los valores de los pixeles representan la reflectancia en superficie en distintas longitudes de onda. Acá pueden encontrar el manual de estas imagenes donde les detallan la descripción tanto de los nombres de lor archivos como de los preprocesamiento que tienen realizados.
Para este ejercicio hemos realizado un clip de cada una de las bandas originales de la imagen y ya multiplicamos a cada una de las bandas por el factor de escala indicado en el manual (0,0001).
Las longitudes de onda y la resoluciones de cada banda de la imágen se describen a continuación:
| Banda | Longitud de onda (nanómetros) | Resolución espacial (metros) |
|---|---|---|
| Banda 1 - Aerosoles | 430 - 450 | 30 |
| Banda 2 - Azul | 450 - 510 | 30 |
| Banda 3 - Verde | 530 - 590 | 30 |
| Banda 4 - Rojo | 640 - 670 | 30 |
| Banda 5 - Infrarrojo cercano | 850 - 880 | 30 |
| Banda 5 - Infrarrojo medio 1 | 1570 - 1650 | 30 |
| Banda 7 - Infrarrojo medio 2 | 2110 - 2290 | 30 |
En la carpeta clip encontrarán los datos que vamos a usar en los ejercicios. Cada banda del clip se encuentra en formato .npy
Importamos las bibliotecas y abrimos algunas bandas de datos
import numpy as np
import matplotlib.pyplot as plt
from google.colab import drive
from mpl_toolkits import axes_grid1
def add_colorbar(im, aspect=20, pad_fraction=0.5, **kwargs):
"""Add a vertical color bar to an image plot."""
divider = axes_grid1.make_axes_locatable(im.axes)
width = axes_grid1.axes_size.AxesY(im.axes, aspect=1./aspect)
pad = axes_grid1.axes_size.Fraction(pad_fraction, width)
current_ax = plt.gca()
cax = divider.append_axes("right", size=width, pad=pad)
plt.sca(current_ax)
return im.axes.figure.colorbar(im, cax=cax, **kwargs)
drive.mount('/content/gdrive')
Drive already mounted at /content/gdrive; to attempt to forcibly remount, call drive.mount("/content/gdrive", force_remount=True).
#Levanto los datos RGB + Infrarrojo
bNIR=np.load('/content/gdrive/MyDrive/Otro/teledeteccion/clip/LC08_L1TP_225084_20180213_20180222_01_T1_sr_band5_clip.npy')
bR=np.load('/content/gdrive/MyDrive/Otro/teledeteccion/clip/LC08_L1TP_225084_20180213_20180222_01_T1_sr_band4_clip.npy')
bG=np.load('/content/gdrive/MyDrive/Otro/teledeteccion/clip/LC08_L1TP_225084_20180213_20180222_01_T1_sr_band3_clip.npy')
bB=np.load('/content/gdrive/MyDrive/Otro/teledeteccion/clip/LC08_L1TP_225084_20180213_20180222_01_T1_sr_band2_clip.npy')
Ploteamos banda roja
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(bR)
add_colorbar(im)
<matplotlib.colorbar.Colorbar at 0x7f164c6fa1d0>
Conviene setear el rango de valores para que se vea mejor. Para saber qué rango, miramos el histograma
h=plt.hist(bR.flatten(),bins=250,range=(0,3))
Usemos entonces el rango [0, 2] para visualizar la banda roja
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(bR,vmin=0.3, vmax=2)
add_colorbar(im)
<matplotlib.colorbar.Colorbar at 0x7fd591ee9b70>
Repitamos el procedimiento con NIR (near infrared)
Por un lado podemos ver muchos pixeles con bajo valor esto es el Río de la Plata, posteriormente una distribución gaussiana correspondiente a la tierra.
h=plt.hist(bNIR.flatten(),bins=250,range=(0,5))
Uso ese mismo rango apra visualizar la imágen.
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(bNIR,vmin=1, vmax=5)
add_colorbar(im)
<matplotlib.colorbar.Colorbar at 0x7f1646ed8090>
Veamos las dos bandas, una bajo la otra
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(bR,vmin=0, vmax=3)
add_colorbar(im)
plt.show()
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(bNIR,vmin=1, vmax=5)
add_colorbar(im)
plt.show()
En la banda del rojo la vegetacion refleja poco mientras que en el infrarrojo refleja mucho.
rgb = np.dstack((bR,bG,bB))
plt.figure(figsize=(7,7), dpi=150)
plt.imshow(rgb)
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
<matplotlib.image.AxesImage at 0x7f1646c97e10>
Si modifico las matrices obtengo imagenes con colores "falsos"
plt.figure(figsize=(7,7), dpi=150)
rgb = np.dstack((bNIR/4,bR/1.5,bG/1.5))
plt.imshow(rgb)
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
<matplotlib.image.AxesImage at 0x7fd718de56d8>
El Índice de Vegetación de Diferencia Normalizada, también conocido como NDVI por sus siglas en inglés. Este índice, basado en la intensidad de la radiación de dos bandas del espectro electromagnético que interactúan particularmente con la vegetación, aporta información sobre la cantidad, estado y desarrollo de la misma. Para calcular el NDVI se utilizan las bandas espectrales Roja e Infrarroja
#calculo NDVI y trunco potenciales problemas
ndvi = (bNIR - bR) / (bNIR+bR)
#el nvdi deberia estar acotado entre -1 y 1, por lo cual normalizamos para evitar errores.
ndvi[ndvi>1]=1
ndvi[ndvi<-1]=-1
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in true_divide
Graficamos el ndvi:
plt.figure(figsize=(7,7), dpi=150)
im = plt.imshow(ndvi,vmin=0, vmax=0.8)
add_colorbar(im)
plt.show()
clases_ndvi = []
for i in ndvi.flatten():
if i < 0:
clases_ndvi.append("No vegetada")
elif (i >= 0) & (i <0.1):
clases_ndvi.append("Área desnuda")
elif (i >= 0.1) & (i<0.25):
clases_ndvi.append("Vegetación baja")
elif (i >= 0.25) & (i<0.4):
clases_ndvi.append("Vegetación moderada")
else:
clases_ndvi.append("Vegetación densa")
Distribución de clases según vegetación
plt.figure(figsize=(10,5))
h = plt.hist(clases_ndvi, bins=np.arange(6)-0.5)
Clases en el mapa:
clases_ndvi = []
for i in ndvi.flatten():
if i < 0:
clases_ndvi.append(0)
elif (i >= 0) & (i <0.1):
clases_ndvi.append(1)
elif (i >= 0.1) & (i<0.25):
clases_ndvi.append(2)
elif (i >= 0.25) & (i<0.4):
clases_ndvi.append(3)
else:
clases_ndvi.append(4)
clases_ndvi = np.array(clases_ndvi).reshape(ndvi.shape)
from matplotlib import colors
import matplotlib.patches as mpatches
# Creo colores
cmap = colors.ListedColormap(['black', 'y',
'yellowgreen', 'green', 'darkgreen'])
# Defino los limites de cada color
limites = [0, 1, 2, 3, 4]
norm = colors.BoundaryNorm(limites, cmap.N)
# Genero leyenda y grafico con leyenda
texts = ['Sin vegetacion', 'Area desnuda', 'Vegetacion baja',
'Vegetacion moderada', 'Vegetacion densa']
patches = [mpatches.Patch(color=cmap(i), label="{:s}".format(texts[i]) ) for i in range(len(texts))]
plt.figure(figsize=(15,8))
plt.legend(handles=patches, bbox_to_anchor=(0.2,1.3), loc='center', ncol=1 )
plt.imshow(clases_ndvi, cmap=cmap, norm=norm)
plt.show()
#elijo un clip de 600x600 pequeño para hacer las pruebas
sy = 1100
sx = 2000
dx = 600
dy = 600
clip = ndvi[sy:sy+dy, sx:sx+dx]
plt.figure(figsize=(7,7), dpi=150)
im =plt.imshow(clip,vmin=0,vmax=.7)
add_colorbar(im)
plt.show()
rgb = np.dstack((bR[sy:sy+dy, sx:sx+dx]/2,bG[sy:sy+dy, sx:sx+dx]/2,bB[sy:sy+dy, sx:sx+dx]/2))
plt.figure(figsize=(7,7), dpi=150)
plt.imshow(rgb)
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
<matplotlib.image.AxesImage at 0x7f162db3d250>
Usamos un algoritmo de clustering para obtener clusters de pixels similares a partir de los datos de NDVI calculados antes.
#importo el clasificador y defino una instancia para clasificar con dos etiquetas
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=2)
#le saco la estructura bidimensional a la matriz NDVI y la llamo datos
datos = clip.reshape(-1,1)
datos.shape
(360000, 1)
#entreno o ajusto el el clasificador con los datos
kmeans.fit(datos)
#usa el modelo ajustado para poner etiquetas
etiquetas = kmeans.predict(clip.reshape(-1,1))
Visualizamos los resultados
Distingue pixeles con y sin agua.
#visualizo los resultado,s recuperando la estructura de matriz
plt.figure(figsize=(7,7), dpi=200)
im = plt.imshow(etiquetas.reshape(clip.shape),cmap='jet')
#add_colorbar(im)
plt.show()
Ahora con tres clases
Tenemos agua, y dentro de la tierra tenemos dos clases: mucha vegetación y poca vegetación
kmeans = KMeans(n_clusters=3)
kmeans.fit(datos)
etiquetas = kmeans.predict(clip.reshape(-1,1))
#visualizo los resultado,s recuperando la estructura de matriz
plt.figure(figsize=(7,7), dpi=200)
im = plt.imshow(etiquetas.reshape(clip.shape),cmap='jet')
#add_colorbar(im)
plt.show()
Con mas clases podemos empezar a ver diferentes grados de urbanización
kmeans = KMeans(n_clusters=5)
kmeans.fit(datos)
etiquetas = kmeans.predict(clip.reshape(-1,1))
#visualizo los resultado,s recuperando la estructura de matriz
plt.figure(figsize=(7,7), dpi=200)
im = plt.imshow(etiquetas.reshape(clip.shape),cmap='jet')
#add_colorbar(im)
plt.show()